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ABdTiUCT 


The F okker-P lane k equation governing the interaction of 
charged particles with fully ionized plasmas, for a certain 
useful energy range of the charged particles, becomes a iirst 
order partial differential equation in space and tine. However, 
the deflection of the par tie.", as is represented by a higher order 
t era . 

It is seen that if the deflection term is neglected, the 
remaining equation can be solved exactly by the method of charac- 
teristics. The effect of the deflection term is then included 
ar a perturbation. To this end, a first-order perturbation 
method is developed in this work. 

It is found that the effect of deflections ic generally 
to increase the energy deposition near the source. In the plasma 
temperature range of 1 to 10 keV, and particle energy range of 
1 to 5 MeV, most of the particle energy is deposited to electrons, 
and the effect of the deflection term on energy deposition rates 
is small. Hence the first-order perturbation method is quite 
satisfactory. However, the deflection term, becomes more important 
as the particle slows down. Thus towards the end of the particle 
range, both the governing equation as well as the perturbation 
method are poor approximations. 



CHAPTER 1 


IHTRODTJCTICH 


3. . 1 F or ew ord 


The manner in which a Tact charged particle distributes 
its energy in a plasma is important in investigations of the 
development and structure of tb oxmonuclear reaction waves in 
."fixtures of deuterium and tritium. 

In an inertial confinement of plasma, fast charged 
particles can he used as an external beam for heating and com- 
pression, or can be produced as reaction products which carry 
a significant fraction of t o reaction energy. The examples for 
such reaction products (alphas, protons, tritons) axe provided 
iy the following reactions fll , 

D+T - (3.5 MeV) + 1 n (14.1 MeV) (1.1) 

D+D T (1 MeV) + p (3 MeV) (1.2) 

^-n + 6 Li ”■ 4 He (2.1 iieV) + T (2.7 MeV) (1.3) 

r;. e re, 1 n is a neutron, and D and T are deuterium and tritium 
nuclei respectively, and p stands for a proton. 

That is why it is important to accurately simulate how 
and where the charged particles deposit their energy. 
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1*2 Background and Literatu re R evi ew 
1 . 2.1 General 

Past ions introduced into a fully ionized plasma, 
either from the outside or as a result of nuclear reactions, 
undergo two distinct processes during thermalization. First, 
they lose energy to the plasma and eventually diffuse in energy 
space until a Maxwellian distribution at the plasma temperature 
is attained. Secondly, they are deflected by the plasma and 
eventually diffuse into an isotropic velocity distribution. 

'the energy loss problem has been studied extensively since it 
is crucial to the energetics ox any plasma system. However, 
tile deflection of the ions has usually been ignored and a 
straight line trajectory during the entire slowing down process 
has been assumed. • 

Several methods have been used to study the charged 
particle slowing down problem. Cooper and Evans [2] have used 
the range-energy loss relation to obtain, by numerical integra- 
tion, the fraction of energy leaking out of a spherical system. 
Antal and lee [3] have derived equations for the energy and mass 
transport of energetic ions using the conservation laws of 
density and energy in phase space. They have replaced the large 
number of small angLe Coulomb scattering collisions by a net 
frictional force. The resulting transport equation resembles 
the one encountered in neutron transport problems and can be 
solved by well— known methods, such as method used by these 
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authors. C orman et al. [4] have numerically solved the particle 
transport equation (F okker-Pl ancle equation) in the multigroup 
diffusion approximation. And, recently, Haldy and Ligou [5] have 
derived a moment method for the. transport of energetic ions in 
aa infinite homogeneous hot plasma, in the absence oi electro- 
magnetic fields and solved the 1 okker-Pl anck equation for plane 
and point sources. 

1.2.2 The F okker-Pl anck Eq uation 

The F okker-Planck equation is mainly a Boltzmann equation 
i till collision term evaluated for small angle scattering [6], 
and is given as follows : 


af 
3 t 


+ 


v- 


V f + 
r 



v f 


v 


( -|-f) 

collision 


(1.4) 


where, £ (r, v, t) is the particle distribution function, and . 

r, v, t) dr dv is the expected 
number of particles at time t, in a volume element dr located 
;.',t r, whose velocities lie in dv about the velocity v, and F is 
the external force, m is the macs of the individual particles. 

In the Equation (1.4) collision term is calculated by 
liosenbluth [7] and expressed with two terms : 


a - Dynamical friction, 
b - Diffusion in velocity space. 
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A useful representation of Equation (1.4) is ^i?en by 
Cooper f8], reflecting the diffusion in velocity snac c, and 
also assuming that the target particle distribution function is 
ucOiwellian. Kidder and Ree [9] have shown that this is a 
reasonable approximation for slowing down of fast ions. The 
P oLker-Planck equation with these approximations can be written 
as . 


"ft + v , m t) 


= S : (r , v, (1 , t) 
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l u < f + st 
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r m o -p 

‘ + — § q(y ) —1— (1-u^) -4-1 

2v m a' a >i 1 ^ ' a p 


(1.5) 


where 


|i 

i(r , v, |i, t) 


S (r,v, |i, t) 


Cosine of the scattering angle. 

The number of particles at time t per 
unit volume in the phase space defined 
by r, v , (i, 

The source term defined in the same way. 


a 

Z 

m 


The plasma component (electron-ion). 
Atomic number. 

Particle mass. 

Mass of the plasma component a, 
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Temp erature, 
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Haldy and Ligou [5] take the energy as basic variable 
P 

(E = mv /2) and introduce flux co as follows : 


cp (r > E, |i, t) = v f(r,E,l-'.,t) 


(1.6) 


where : 

f(r,E,p,t) = f(r,v,p,t) v 2 ~ 


or alternatively it can be written as ; 

„2 

<p(*,E,p,t) = — f (r , v, p,, t) (1.7) 

l’hen for plane geometry, the Equation (1.5) takes the form : 
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For slowing down domain, the following additional assumptions 
■-an be made [5] ; 

1. E >> 0 ± 

2. E << ~ e 

in e 


('these conditions are fulfilled in a large domain, e. 3*5 MeV 
alphas in 5 KeV plasma, 1 MeV proton in 5 KeV plasma). With 
these assumptions, Equation (1.8) takes the following form s 
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end Coulomb logarithms can be calculated from Ann end ix I 
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1.3 Outline of the Present Work 


In Chapter 2 some mathematical methods which are utilized 
in this work for the analytical treatment of the approximate 
1 okker-Planck equation are reviewed. Also the zeroth order solu- 
tion of the P okker-Planck equation has keen given in that chapter. 

The first order perturb. tion method is then presented in 
Chapter 3. In Chapter 4, the nu'erical results obtained are 
resented and discussed in graphical and tabular form. Nummary 
and recommendations for future vork are included in Chapter 5- 
Come detailed calculations and properties of some special function 
are given in Appendices I to I^. 



CHIP TER 2 


MATHEMATICAL BACAORC/UHD AND THE ZEROTxi OEDLL SOLUTION 

2 . 1 The Method of Characteristic s 

The E okker-Planck equation, Equation (1.9) s is a linear, 
second order partial difi er enriai equation in four independent 
variable (X, E, p., t) . But one e straight line traj ectory has 
beer assumed, the equation will be reduced to a n on-homogeneous, 
linear, first order partial differential equation in two indepen- 
dent variables (X and E) for tine independent cases. Thus, 
Donation (1.9) will take the following- form , 


3cp(X,E, p, t) 

3 X 


(0(B) d ) 
3E 


S(X,E,p,t) 


( 2 . 1 ) 


Since this form of the equation will be used very frenuently 
during, both the zeroth order and the first order solution of the 
I oJ::-.er-Planck equation, we will outline the general solution of 
t.-e non-homogeneous, linear, first order partial diii cr ential 
eolation in the following linos. 

In general, a non-homogereous, linear, first order partial 
differential equation can be written in the following form : 


a(x,y) u x + b(x,y) u y 


c(x,y) 


( 2 . 2 ) 
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wh er e , 


u 


x 





9,u 

9 y 


a, Id and c are coefficients and they are the functions 
of x a.nd y. 


If 

u = u(x,y) 

is a solution of Eq. (2.2), at a point x, y, u of su? face 
u(:t,y), the vector (u , u , -1) has the direction of the normal 

x y 

to the surface [12] and Eq. (2.2) states that the vector (a,b,c) 
is perpendicular to the vector (u^, u , -1) and must lie in 
tangent plane of the surface u(x,y). That is why any solution 
surface u(x,y) through the -joint (x,y) must Le tanpent there to 
prescribed vector, namely, to the me with the components 
a(x,y), b(x,y) and c(x,y). So, instead of solving a partial 
differential equation, if one finds a surface u(x,y) which con- 
sists of curves tangent to vector (a,b,c), then u(x,y) will be 
the solution of Eq. (2.2). The curves in the surface u(x,y) 
arc called characteristic curves and should satisfy the following 
relation ; 


dx 

a 


b 


du 

c 


(2.3) 
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(Proof can be found in [13]). In this way a linear, iirst 
order partial differential en nation can be converted to two 
ordinary differential equation defined by Equation (2.3) which 
can eas'ily be solved. 

2.2 The P erturbation Method 


The perturbation method is an approximate method which 
is commonly used in many fields of science when enact solutions 
oi the equations are not possible. In this method it is assumed 
i at the dependent variable cm be written as the sun of two 
,-:r ts . One of these parts is ox sufficiently simple structure 
which can be easily solved, and the other part is small enough 
o that it can be regarded as a n erturbation on the urst part. 

This method can be very well applied to the Fokker- 
planck equation to study the effect of the deflection term on 
the energy deposition. For this purpose Eq. .(1*9) can be written 


3Uco 

3X 


U 


iEl _ ! i- (1-b) 112 

ah 1 3|-l V • J 3 P- 


U S 


( 2 . 4 ) 


or 


with $ 


U cp 




3$ 
3 X 


T1 1±_ T i_ n-u 2 ) 

U 3E T 3p ^ ] 


U S 


(2.5) 



12 


i.ere we define the following two operators, H Q and E-^ . 

H 0 = ^-UCEjr- (2.6) 

H 1 =. -fn d-P- 2 ) (2.7) 

In terras of these operators Equation (2.5) clik be writ'ten as: 
(h Q + %) *= US (-.8) 

ft. apply the perturbation met: ou we can writ-, n two parts : 

* = ® 0 + (2,9) 

where is a sm all perturbation on $ . 

If Equation (2.8) is solved as it is, it will give the 

itiul $ , but this is not p ossicle at this stage. If it is 

civ . d by neglecting the deflection term (E^ $ = o) it will 

• ive *. which we call the reroth order solution. If it is 
0 

solved by replacing * by $ 0 + solution will give the 

-j (by neglecting the second order terms) which we call as the 
: irst order solution. So the seroth order solution will be the 
solution of the equation : 
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H 0 $ 0 = U S (2.10) 

whereas the first order solution will be that of the eolation, 

11 o 5 1 - - H l ? 0 O- 11 ) 

In the following section the solution of the Donation 
(2.10) and in Chanter 3? the solution of the Equation (2.11) 
will be given. 

2.3 The Zeroth Order Solution 

The solution of the steady state Fokker -Planch equation 
: j neglecting the deflection term (T(E) = 0) will be named as 
the zeroth order solution. It is nothing but assuming e, straight 
line trajectory during the entire slowing down process. With 
this assumption, the governing cruation will take the following 
i onu . 

ct 4 _ ^ $ 

_2_ „ U(E ) _o_. = U(E) S(X,E,n) (2.12) 


vn ere. 


© 0 = $ 0 (X,E,|a) = U(E') <p 0 (X,E,n) 

9 c (X,E, t-0 - Zeroth order angular flux. 
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The solution of the above equation is easy by the method of 
ciiaracteristics, details of which are explained in Section (2.1). 
!■' or this case the characteristic equations will be as follows : 


dx _ dE 
-U ( E) 


d (Z , E, p.) 

U(E) s (z ,e,V) 


(2.13) 


first part of the above equation gives 

= X(E) - X (E* ) 


where , 


X (E) is defined by ; 


E 


s 


X (E) 


/ dE 
U(B) 
E 


(2.14) 


(2.15) 


where E^ is the maximum energy of the particle, physically x (E) 
is the distance travelled by a charged particle while its energy 
is reduced from E 0 to E by continuous energy loss [14]. This can 
easily be seen from the definition of U(E) which is the total 
continuous energy loss per uni"., path length from small angle 
scattering. U(e) has a unit of energy per distance, whereas X (E) 
is in distance units. 

The second part of Eqs. (2.13) can only be solved for a 
specified source tern. Throughout the work a plane mono energetic 
isotropic source will be used in our calculations. This source 
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coi'i be written as 2 


S ( :c , E, p) = 6(x) 6(E g -B) 


( 2 . 16 ) 


ir,iere : 


6 is the Dirac delta function, and 
S Q is the number of particles produced per unit 
area per sec. 

how the second part of the Eqs. (2.13) can be solved with this 
s un rc e 1 

B 

$ 0 = - / S(x>, e’, p) dB’ (2.17) 

B s 


With the source term as in Eq. (2.16), we obtain 


0 

o 



/ 6 (x ' ) 6 (Eg-E ' ) dE' 

E 


(2.18) 


•’here x 1 is defined with Eq. (2.14), so that it should be 
replaced in terms of e' : 

S s 

Cj 

= j | / x(E) -!- p. X (E ’ ) ) 6 (Eg-E* ) dE ' 

E 


(2.19) 
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The above integration can be solved, by using the properties 
01 the Dirac delta function (see appendix II) : 


S 


o 

4tc 


6(x - p x (jj) + p x(E') ) i 


( 2 . 20 ) 


;E =E, 


_ o 
4tc 


6(x - p ' x (E) ) 


( 2 . 21 ) 


iron solution, the angular flux density can directly 

be written ; 


cn Q (X,E,P) 


U(B) 


o 


4itU(E) 


6(x - p x(E)) 


( 2 . 22 ) 


And flux density can then be calculated from the zeroth order 
angular flux density by integrating over p; 

1 

to (x > B) = 2“ / rsfii 6(x ° ) itl (2 - 25) 


for x > 


o 
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wh er e : 


x Q = x - p x(E) 

dx Q = - x(E) dp, (x and E constant). 


I 1 !. '.us Equation (2.23) can be rewritten as ; 


x- x ( e) 


cp n (x,E) = - 


S 


- f 6(x ) 

2 U(E) X (E) 1 ° 0 


x+ X (e) 


(2. 24) 


9 (x,E) = 


S 


x+ x(E) 

/ 6 (x ) dx (2.25) 


2 U(E) X (E) ' ' 0 0 

x- x(E) 


S 


2 U(E) X (E) 


if — <_ 1 (2.26) 


x| 


0 


if 


7CI7 


> 1 (2.27) 


finally the expression for the rate of energy deposition can be 


written as : 

E S (x) 

(x) = / U(E) n 0 (x,E ) dE 

B mi*i 


( 2 . 28 ) 
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Using cp Q (x,E) from Eqs. (2.26) and (2.27), we get i 

E*(x) 

/ 1° dE (2,29) 

2 X (E ) 

E . 
mm 

where E S (x) is the maximum energy of the particle which it can 
have at a position x. This is the energy, from the definition 
oi x ( E) , where 

•X(E* ) = x 

In the other hand, E min is the minimum energy of the particle 
nhich it can have at a position x. Although it seems that 
should be. zero which is the minimum energy, first, it can not be 
less than the plasma temp era tune since particle will not lose 
energy below the plasma temperature; secondly it can not be less 
than a specified value of energy above plasma t emp era cur e since 
the Equation (2.1) has been derived under the assumption indicated 
bf Equation (1.8a) and it will not be valid below a specified 
vrlue. This energy value has been taken as ten times the plasma 
temperature through most of the present work. The energy deposi- 
tion can be written separately for ions and electrons ■; 

B*(x) 

ifpx) - / -|- f, o (x ’ B) dB 

"min 



(2.30) 
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(2.31) 


( 2 . 32 ) 


(2.33) 

It can be easily seen that the maximum range, x , of the 

it. 

particles will be 

\ = *<W’ i 10 e - 

Thus Eqs. (2.28) to (2.33) are valid ± or x < x^. For x > x^. , 

energy depositions are zero. 

In the next chapter the first order solution of the 
1" okher-Planck equation will be derived by applying perturbation 
uethod as approximation and the method of characteristics will 
bo utilized for the solution of the differential equation. Nume- 
rical results based on either the zeroth order solution or the 
first order solution, will be presented and compared in Chapter 4. 


I 1 (x) 
PI 


h (x) 


/ 


b S 


2 u("e) x(e7e dB 


E . 

mm 


b e i« 


B S (x) 




a h- co 0 (x,E) dE 


m 


E* (x) 


a E^ s o 


2 U (E) X(E) 


dE 


E . 
rmn 


CHAP mi 3 


THE PERTURBATION METHOD 


3 » 1 The Eirst Order Per tur ha. t i o n Method 

In the preceding chapter it has heen shown that the 
iirst order solution will he the solution of Equation (2.11) 
w ich is : 


H o = - H, o Q 


or, 


3 $ 30 

pr - tf 


T(E) ^d-^ 2 ) * (B) h 


(3.D 


or, 


3»- 




3 X 


U 


3 0- 

aE 


S.T o 2, 

— fl-u 2 ) 

An U ^ ' 2 

3 p 


— 2— •( 2u)**--^- 
An K 3 ii 


(3.2) 


S 1 (s, E, \i) 


where : 


T = T(E) 


U = U(B) 


6 


6 {X - (i x(E) } 



21 


I‘u is easy to solve this equation with the following characteris- 
tic equations : 


dx dE 

v- ~ u (e) 


d $ 


1 


°11 + S 12 


( 3 . 3 ) 


where : 




S o T ^ (1 2, a 2 6 £ x-u > (J]2 

4,1 


S-^2 — ®i2 (x > E> M') 


S o T , v 36 { X ~V x(B) } 

-45T- (2**) — 


and 


’1 


S 11 + S 12* 


From the first of the characteristic equations, Eqs. ( 3 . 3 ), we 

have, 


x - x * 


X(B) - X(B') 


(5.4) 


Iran the second of the Equations (5*5) we have, 
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E 

> S> p) — — j 

E g 


&ll(x ' ,E* , p) + S 12 (x' ,E' , !-0 


U(E') 


dE 


(5,5) 


there s' is feiven by equation ( 5 . 4 ) and should be replaced in 
terms of E f before carrying out the 'dE' integration. For simpli- 
city let. 


'1 


‘11 + *12 


(5. 6) 


where : 


E. 


S. 


'-'ll (■£ » E, 4) 


/ dB - 


E 


E, 


E 


/ ^-5 6 { x'-p x(E)} dE* 

J u ( m ) a p 

(5.7) 


a P 


where, it can easily be shown that the Dirac delta function is 
independent of E 1 by replacing x ' in terms of E ’ , x and p from 
Equation (5.4). So, 


6 {x' - p X (E* ) > = 6 { x - p x (E) } 


(5.8) 


Hence this term can be taken out of dE' integration as a constant 
term. Thus we have, 
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e l;L (x,E, n) = 


S 

4it 


° (1-P 2 ) I 

3 P' 


(3.9) 


W''ere: 


1(E) = / 


E 


, dE 
U(E ) 


Similarly , 


E„ 


<?2_2( x » I 1 ) - / 

E 


Sn o(X , E j (-0 T 

r dE 

U(E ) 


k (2 ^ li 1 


(3.10) 


Erctn Equations (3«9) and (3.10), the angular flux density can 
directly be written as ; 


fi 


?u + <p 12 


with, 


®ll(-^ t E> 1-0 


S rt o A 

2 (,l-|i 2 ) 1A j; 

4tc U (E) 3 ja 


c f 1 12^ :x:,, ® , ^ 


S 


4% U(E) 


<») If 1 


(3.11) 


( 3 . 12 ) 


(3-13) 
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,.>j- performing the dp integration the flux density can be obtained 
as follows ; 


<Ph(x,E) - 2tc j q) (x,E, p) dp 


>jO, 


S 0 1( b) 

2 U(E) A (E) 


- 2 


*7 / (l-l 1 ) 6(x ) dx. 

A Ep^ 0 


(3.14) 


wrere i 


x n = x-p x(E) 


dx = - X (e) dp \ (with x and E constant) 


co ll ( :c>E ) = 


S Q jjg) 3 2 (1-P 2 ) 


2 U ( E) X(E) 3p‘ 


(3*15) 


x o =# 


(See Appendix II for the properties of the Dirac delta function) 
If Equation (3*15) is evaluated the flux density will he . 


9ii( x »®) 


S 0 1(E ) 
U(E) X (E) 


(3*16) 


Similarly it can easily he shown that : 

S o I(E) 

cp 12 (x,E) = - U(E) X ( S ) 


(3.17) 
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Thus the total angular flux density will he 2 


cp(x, E, p) = 0 (x,E,p) + w-i (x, B, (J.) 


(3.18) 


cp (3C, E, M-) = 


rs + (l-sT) 1 ~ 2 - (2M) I ] (3.19) 

3 \X 47 


whereas the total flux density will he the same as the seroth 
order solution ; 


cp(x,E) = cp Q ( x , E ) + cp (x, E) 


( 5 . 20 ) 


co(x,E) = 


S 


+ 


S Q 1(E) 


2 U(E) X (E) U(E) X(E) 


S Q 1(E) 
X(E)U(E) 


o 


x (e) 


(3.21) 


Thus it is seen that the perturbation source affects the 
angular flux density, whereas it does not affect the flux density. 

This is also true for the higher order perturbations. In 
Appendix III we have shown that the second order perturbation will 
lead to the same result. In the next section we will point out the 
reason for ti.is behaviour and accordingly develop a modified first 
oruer perturbation method to see the effect of the deflection term. 
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- • 2 The nature of the Deflection Term 

In the preceding section it has been seen that if -h^ 
term of Eq* (2.11) is used as a source term in the fir at order 
perturbation method, it does not affect the flux density, and 
hence energy deposition rates. In this section we will look at 
the perturbation source term S- L (x, E, p) in more detail. It can 
be written as sum of the two terms ; 

S 1 (x,E,p) = S i;L (x,E,p) + S 12 (x,E,p) (3.22) 

To see the nature of the source it can be- integrated over p : 

1. 

Sn(x,E) — 2k ! 3 1 (x,E,p) dp (3.-23) 

-1 

with, 

x Q = x - p X (E) 

and 

dx 0 = - Me) dp 

it will be : 

s u U.B) = -IfiT / s x (x ’ E ' rt dx o 



27 



This clearly shows that the source of the first order perturbation 
equation, Eq. (2.11), is not introducing any new particles to the 
Ojstem as it is expected. Instead, it is redirecting sane of the 
particles from their original path. And this effect is such that 
at each, x and E, S 1]L /CJ(E), particles per unit time, per unit energy 
interval and per unit volume, experience deflection. This is the 
overall effect of the deflection term, need not he true for a sin- 
pie particle to he deflected 180 degrees. That is why in order to 
ho able to see the effect of tile deflection term, for oh e source 
term S 12 (x,E,li), the E okker-Planck equation should he solved in 
the reversed (— p) direction. The solution for the soi.rce term 
has already been obtained in Section (5.1). In the 
next section the solution ton ft- At m nl 


*h Pol onl 
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2 * . T hejyL ^odi f j ed First Order Perturbation M eth o d 

The same procedure of Section (3.1) will be followed to 
xi-id the modified perturbation solution with the modified source 

t orrn : 


ju S T(E) 

^12 = ~ ( ^h) 'a]7 ^ { X+ P l *-(e) i (3.27) 


defined in the range of x between 0 and x^. And with this source 
term the equation to be solved now is s 


P' 



a x 


- 13(B) 


m 



»E 


q m 

°12 


( 3 . 28 ) 


which can easily be solved with the following, characteristic 
conations : 


dx _ dE 

[a “ ~ "tfOiT 


m 


d C 


12 


m - 

W 12 


(3*29) 


Again, as in Section (3.1), the first of Eqs. (3*29) givess 


x - x 


X (E) - X(E') 


(.3.30) 


where the second gives : 
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m 

$ (x,E,p.) 

12 


E 


1 


S‘ (x,,E ,,p) 

/ — — : dE 1 for x 

U(E ) 


m 


o 


E 


(3.31) 


Bn 


-m / ' * 


t Sip(x , E ,p) 

J — hd :■ dTC 


E 


T J(E ) 


dE for -x x /. o 

ill 


(3.32) 


where the limits can easily be decided by considering the minimum 
•ihd maximum energy of the particle which can have contribution to 
'the flux at any x and E. It should also be kept in mind that to 
lave a particle at energy E at any x, the total distance travelled 
by the particle starting from source energy E 0 should be same, vis. 

o 

X(E). With this consideration E^_ can easily be calculated as : 


\(E 1 ) 


x + x ( e) 


(3.33) 


udng the actual expression for fa® 2 (x,E,p) for x > o ; 

Bn 


m 


0 (x,E f n) 

12 


S o j T ( E * ) - 3 


. . 7^ 2p ~ 6 {x'+p X(e') > dE* 

4lx E U(E )' ^ 


(3-34) 


where x* is given by Equation (3-30). So using this value of x’ 


in the above equation : 
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E 


U,B,p) 

12 


S 


1 


Of 2(E) 3 

4* 1/)' 2i ' 1 ^ 6 {x_d X(E)+2 ^ *(!)> as' 


E 


(3.35) 


int ©grati on can be rjono _ 

e aone '~ nu angular flux density can be 

:tcv ,d. But since our train interest is to find the flux density 

' Cner& dep0siti on rates, „ rtll directly CTt J.n ; .te the 

'O'- integration in order to po~ : ~ -i-u a -f-i n , v a ■ . 

* g Gu Cue flux density since it is simp- 

i.or to carry out the dp integration. The flux density is : 


uith 


m 

* (x,E) 

12 


*1 

I clE ’ 
E 


/ ie’ r in : i 


-i 


X - p x(E) + 2p X(e’) 


C(/( 2 "Ib 6( y o> 

(3.36) 


dy. 


{2 x(E ) - A (33)} dp 


hniu 


m 

< (x,E) 

12 


En 


8 , 


/ d B ' 


m 


E 


...Hi 2 . jl_ 

U(e’) 


6 (yJ 


d 7 


o 


F 2 x(s )- A(B)] 


(3.37) 
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I: sing the properties of the Dirac delta function given in 
Ann end ix III, dp integration can be carried out easily giving 
tl ; e 1 o 11 owi rig res ul t : 


m 

(x,E) 

12 


s , 


/ 


ME 1 ) dE 

v(e'T [2 x(b')- Me)] 


y )■ X > O 
FI " 


(3.38) 


1 • i ,11 arl y i or th ere gi on -x . 

J c in 



g" x £ o , 

•v ’ 7 

T(e') dE ' 

»■ hi iw W m mm rn m m m mm * »■■„■■■ i nn mm mm* m m — 

U(e') [ x(E)-2 X (E )] 

(3-39) 


fans 


the perturbation flux density can now be written a.‘ 


rn 


•m 


^(XjE) = 0 12 {x f 'S) -I- 5 12 (-x,E) - 0 1]l (x,e) 


(3.40) 


In the above equation although - 1 (x,E) is the flux m the posi- 
tive x direction the term o° ? (~x,E) has been added since there 
will be an equal amount of contribution to the flux ox x and E 
■from the distributed source in negative x region as a result of 
f e solution of the zeroth order equation in that region. That 
is why instead of solving the same equation in negative x direc- 
ti on, :“ 2 ( —x , E) has been used considering the spnsetr^^^he 

•■iroMem. 


<**> 


; 

711 ^ 1 ^ ^ 


^ \ 
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The first order flux density is then, given by the folio- 
vin t equation ; 

? (^,E) = d 0 (x,E) + c 1 ( io F B) (3.41) 

or, 

, . ^' 0 (x,E) (x,E) 

cp(x,Bj + (3.42) 

U(E) U(E) 

lie energy deposition rates can be calculated from the Equations 
(2.28), (2.30) and (2.32) simply by replacing the zeroth order 

density with the first order flux density for all the cases* 
These energy deposition rates have been calculated numerically 
end the result will be given in the next chapter, together with 
the results for the zeroth order solution. 



CHAPTER 4 


NUMERIC A! CALCULATION’S AND RESULTS 

4*1 Energy Balance at Steady State 

The total energy deposition rates of zeroth and first 
ora or approximations have beer calculated numerically using the 
analytical equations derived in Chapters 2 and p* In order to 
: sure about the result obtained an energy balance has been 

<. *n e. 

T<. tal energy deposition rate as calculateu from either 
Mi ci zeroth order flux density or from the first order flux den- 
rxty must be equal to the total energy introduced by the source 
nor unit time. The total energy introduced by the source per 
iu-.it time ner unit area for a plane source can be written as ; 

1 

S 

E T = 2n / d|i / E dE / dx ~~ 6(x) S(E q -E) (4.1) 

-1 

= S '$ (4. 2) 

o s 

however, the energy deposition of the particles has been studied 
upto a minimum energy of ten times the plasma temperature. Hen- 
ce, actual total energy deposition rate per unit source area for 

Kali* 3pace will be s 
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■%> ~ 2 (E g - 10 0) (4.3) 

Total energy deposition rate per unit source area iron the 
"eroth order approximation is : 

x 

m 

Ema = f I (x) dx (4.4) 

xu pi 

0 

whereas that of the first ordei approximation is ; 

x 

m 

E T1 = / Ip-L 00 dx • (4.5) 

o 

Equations (4.3)? (4.4) and (4.5) have been calculated numeri- 
cally for various plasma and charged particle properties. The 
results are tabulated in Table 4.1 for the deposition of 1 MeV 
proton energy to a DT plasma at solid density for various plasma 

temperatures. 

Similarly, the same energy balance has been done for 
'•joint source. Instead of finding analytical solution of the 
point source problem, the energy deposition due to a. point sour- 
ce is d irec tly calculated from the plane source to point source 
conversion relation [15] which is : 



S55 
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W r) 


1 

2kt 



j 

j 

\ x=r 


(4. 6) 


Energy introduced from a point source is twice of that of the 

"•lane source given by Equation (4.3), since for point source 

full space will be considered. Energy deposition of the point 

source in zeroth order approximation is ; 

r 

m 

E tq = / I p ° t (r) 4ur 2 dr ( 4 , 7 ) 

o 

where s 

r ffi = maximum ran/ e of the particle in radial 
direction 

= energy deposition of a point source in 
zeroth order approximation. 

Similarly for the first order approximation total energy 
deposition is : 

r 

m 

E I1 = / ^t ( r ) 4ltr2 dr (4*8) 

0 

The quantities E^, E TG and E T1 , for plane and point 
isotropic sources, for three different plasma temperatures 
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al e listed in Table 4*1. It can be seen that both the zeroth 
order and the first order flux densities calculated in tliis 
work satisfy the energy balance requirement to a very high deg- 
ree of accuracy. 

d . 2 Energy Deposition Calculations 

For the calculation of spatial energy deposition rates 
based on the zeroth and first order approximations developed in 
Chapters 2 and 3> a computer program, the details of which are 
given in Appendix IV, has. been written. For various energy 
intervals, the accuracy of the numerical integration, differen- 
tiation and energy balance has beon checked. As a result 200 
energy intervals and simple rectangular rule has been used dur- 
ing calculations. 

The computer program has been run for different plasma 
temperatures, changing from 1 KeV to 50 KeV; for two different 
olasma densities, solid density and thousand times of the solid 
iensity; and for different charged particles with different 
energies, mainly 3.5 MeV alphas and 1-2-5 MeV protons. Results 
of those calculations will be presented in the following sub- 


sections 
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4 - 2 ‘ 1 Deposition in Zero th Order Approximation 


| Tlle energy deposition in zeroth order approximation is 

: presented in graphical form in Pip. 4.1 and Pig. 4.2 for 1 MeV 

I ■’ 1 ‘ 0t0n8 and 3 ’ 5MeV ^P^s, respectively. A dimensionless q uan- 
; tity E^U) x m /S o B g is plotted versus a dimensionless quantity 
| --/my Each graph includes the Total energy deposition as well 
uhe energy deposition to ions and electrons separately. The 
.nun characteristic of these g.aph.s is that the energy deposition 
to electrons is much more than that to ions. This feature has 
■,eeri observed for the first order solution also; and lor both 
'lune and point sources. To s o this feature clearly, results 
are also t iv en in tabular form in Tables 4.2 and 4.3. 

To find the reason for this behaviour, which in our 
hn owl edge has not been mentioned in the available literature, 
vc studied the effect of different plasma temperatures, plasma 
densities and charged particle energy. As a result we have seen 
u ' kt mainly the plasma temperature is affecting the oenosition 
1. 0 ions and electrons. This can be seen analytically a,l so. 

Energy depositions to ions and to electrons are propor- 
tional to b/E and a) E respectively, since interactions with ions 
1m represented by b/B and that oi electrons by aYfi. The ratio 
of the energy deposition to ions, ^ (x), to the energy deposi- 
tion to electrons, ^(x) will be proportional to b/a for a 
specified energy which can easily be derived from Equations (1.9a) 
..uci (1.9b). The plasma temperature strongly affects the ratio 




X / X 

Pt9* 4'1 Deposition of 1 MeV nrntnnQ 




T 6&rp er a Pup # ( f op pi ap. o 
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b/a sine© both b and a are the functions of the plasma tempera- 
tui e directly or through the’Debye length. Table 4.4 gives the 
v.‘li, es of a and b for various plasma temperatures. As a result 
b can be concluded that for low plasma temperature (1 to 10 Sev) 
u .. e energy deposition is mostly to electrons whereas at higher 
u«..'iperature (50 K'ev) the energy deposition to electrons and ions 
uill be comparable. 

'' ' 4 * ^ Energy Deposition in the first Order Appr oxiaa ti on 

In figures 4-3» 4.4 and 4. 5, the solution of the first 
ora or approximation is plotted together with the zeroth order 
solution for total energy deposition as well as the energy depo- 
sition ions and electrons for a plane source of 1 1-ieV protons 
ixi 10 kev Df plasma. Again as coordinates, the dimensionless 
m.c'in titles are used. Comparison of the zeroth and the first 
order solution shows that both (solutions are very close to each 
other. That is what is expected since at 10 Kev plasma, tempera- 
cure most of the energy is deposited, to the electrons rather 
chon ions which can be seen frcm the comparison of the figs. 4.4 
■aid 4.5. Since the deflection oi the charged particles is mainly 
i£ cause of the ions, for low t cm eratur e plasmas, deflections 

iave negligible effect on combined deposition to ions and elec- 

;rons. 

In Fig. 4.6, the energy deposition of 3.5 heV alphas 
o 50 heV DT plasma shell of radius r for a point source is 
lotted using dimensionless quantities. Plasma density is taken 
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TABLE 4.4 


Effect of the Plasma Temperature on the Coefficients 
a and b in Equations (1.9a) and (1.9b). 


Plasma Temperature 
[ KeV ] 

a 

[ (J )Vm] 

b 

[(J) a /m] 

1 

0. 473xlO~ 2 

0 . 203 x 10" 24 

5 

0. 580x10" 3 

0. 288x10" 24 

10 

0. 2 29x10" 3 

0. 3 24x10" 24 

50 

0. 246x10" 4 

0.490x10" 24 





>f DT plasm 





Plasma at 50 KeV and 1000 Solid Density 



4 6 Deposition of 3-5 MeV alpha energy to DT pi 
shell of radius x for a point source 
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;is 1000 times the solid density in order to he able to compare 
our results with the available literature. The zeroth order 
solution is the solution of the reference [51. Our solution is 
:.'cre meaningful since at the shell of radius r , the energy 
icposition is expected to be zero, not infinity. Also this first 
order energy deposition curve is more similar to the energy depo- 
sition curve of Haldy and Ligou [5] in which the approximate 
tolution of the Fokker-Planck equation has been obtained by a 
moment method and the solution has been found upto 11 ^ Legendre 
polynomial coefficient. For comparison the solution oi the 
moment method and the first and zeroth order approximate methods 
is plotted in Fig. 4-7 for 5.5 heV alpha particles, in 50 KeV and 
thousand times the solid density D1 plasma for a. point source. 

In Fig. 4.8, the energy deposition to ions and electrons are com- 
pared separately. It is seen that the energy deposition to 
electrons shows a very good agreement in both cases whereas for 

r h G. 6 x there is a considerable difference in energy deposi- 
m 

tion to ions, because of the deflection term. 

4, 2. 3 Effect of Plasma Temperature and Density 

In Fig. 4.9, energy deposition of 1 MeV protons in a 
solid density DT plasma is plotted for two different plasma tem- 
peratures. This plot can give an idea of temperature effect on 
the uotal energy deposition. However its, main effect is on the 
relative energy deposition to ions and electrons, as we mentioned 


earlier. 
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In order to see the effect ol the plasma density the 

solutions for 1 MeV protons in 10 I;.eV D r i plasma are f ound for 

two different plasma densities (solid and thousand times of tae 

solid density). The results obtained are presented in Tables 

(4.5) and (4.6). It can he seen that the dimensionless energy 

deposition, L,-, (x) x /S n E 0 , versus the dimensionless distance, 

r:/x , is virtually the same for the two densities, tote in zeroth 
' in 

order and first order approximations. 

In Chapter 5> a summary of this work is presented 
cstating main conclusions. 
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TABLE 4.5 

Effect of Plasma Density on Dimensionless Energy 
Deposition in First Order Solution 

(Energy deposition of 1 M eV proton in 10 KeV DT 
plasma for a plane source) 


X /x 

7 m 

Solid Density 

1000 x Solid Density 

0.993 

0.138 x 10 -2 

0.141 x 10 -2 

0.855 

0.411 x 10" 1 

0.414 x 10 -1 

0.727 

0.958 x 10 _1 

0.965 x 10" 1 

0.455 

0. 294 

0. 295 

0. 230 

0.651 

0, 652 

0.033 

1.832 

1.836 



TABLE 4* 6 


Effect of plasma Density on Dimensionless Energy 
Deposition in Zeroth Order Solution. 

(Energy deposition of 1 M eV proton in 10 KeV DT 
plasma for a plane source) 


x/x m 

333 

Solid Density 

1000 x Solid Density 

0.993 

0.2257 x 10“ 2 

0.2257 x 10 -2 

0.855 

0.5107 x 10* 1 

0.5111 x 10” 1 

0.727 

0.1081 

0.1083 

0.455 

0 . 3035 

0. 3041 

0. 230 

0. 6478 

0. 6496 

0.033 

1.7795 

1.7851 



CHAPTER 5 


SUMMARY AND CONCLUSIONS 


5 • 1 Summary of the Present "Fo rk 

In Chapter 1 the importance of the energy deposition 
of energetic charged particle?- is explained. The recent lite- 
rature is reviewed together with the F Okie er -Planck equation. 

The general as well as the approximate form of the F okker-Planck 
equation with the main assume tions is presented. 

If in the F okker-Planck equation the deflection term 
is neglected, the equation will be in the form of a linear, 
first order partial differential equation. The general method 
of solving such an equation and the solution of the approximate 
T okker-Planck equation with an isotropic, monoeneip: ctic plane 
source is given in Chapter 1. 

In Chapter 3, instead of neglecting the deflection 
term, a perturbation method is used to solve the 1' okker-Planck 
equation. The nature of the perturbation source is discussed 
and the first-order perturbation solution is obtained using 
the method of characteristics: 


cp 0 (x,E) 


2 U(E) X (E) 


(5*1) 


=■ cp 0 (x,E) + 0 1 (x,E) /IKE) 


cp(x,E) 


(5.2) 
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where ^ is given by equation ( 3 . 40) and U(E) and X (E) are 
as in Eqs. (1.9) and (2.15) respectively. S 0 in Eq* (5.1) is 
the number of particles per unit area per second. 

Once the flux densities are known, the total energy 
deposition rates, as well as energy depositions to ions and 
electrons, separately, are easily obtained. 

Numerical calculations based on the analytical expre- 
ssions developed in Chapter 2 and Chapter 3 are then presented 
in Chapter 4. 

5.2 Conclusions and Recommendations 

The following conclusions can be drawn from the 
present study i 

(a) The effect of deflections of the charged 'particles 
as they slow down in a plasma is generally to increase the 
energy deposit on rates near the source, and accordingly, to 
decrease the same towards the end of the particle range. 

(b) For the plasma temperature range of 1 to 10 ReV and 
the heavy ion energy range of 1 to 5 MeV (roughly this is the 
range over which the F okker— Planck equation used in this work 
applies) , it is seen that the effect of the deflection term 
is small (within 10 % or so) for 

0 < x < 0.6 
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where x m is the maximum ran.;. e oi the charged particles, 
hence it can be concluded that in this range, the first-order 
perturbation method developed in this work is quite satisfac- 
tory. 

(c) The effect of the deflection term is less important for 
the charged particles of higher energy (e.g. for 5 k eV protons 
compared to 1 MeV protons) , and for heavier charged particles 
(such as alphas compared to protons). 

As is clear from the above conclusions, the first-order 
perturbation method developed in this work does not give good 
results towards the end of the particle range. In fact in this 
range, the approximate form of the F okker -Planck equation used 
here is no longer valid. It will, therefore, bo necessary to 
develop more suitable techniaues based on an appropriate govern- 
ing equation. 


5 

I 



REFERENCES 


[1] S. Closstone and R.H. lovberg, "Controlled Thermonuclear 
Reactions", Van Nos tr and , I960. 

[2] R. Cooper and E. Evans, Phys. Fluids 18, 332 (1375). 

[5] M.J. Antal and C.E. Lee, J. Ccrnput. Rhys. _20, 29S (1976). 

[4] E. G. C orman, W.E* Lowe, G. Cooper and A. Winslow, Fuel. 
Fusion, 15, 577 (1975), 

[5] P • A. Haldy, J. Ligou, Fuel. Fusion, 17.6, 1225 (1977). 

[61 J.L. Shohet, "The Plasma State", Academic Press, 1971. 

[7] M. Rosenbluth, W. MacDonald, and D. Judd, Piiys. Rev. 

107, 1 (1957). 

[C] G. Cooper, Lawrence Livermore Laboratory, California 

UC ID-16157 (1970). 

[SI F. Ree and R. Kidder, An analytic solution of the 

therm alization of a fast ion in a plasma, TJCR1-7025 
(1962) . 

[10] N. A. Krall, A.W. Triv clpiec e, 'principles of plasma 
Physics", McGraw-Hill, 1975. 

[11] L. Spitzer, "Physics of Fully Ionized Gases", Inter- 
science, 1962. 

[12] P.R. Garabedian, "Partial Differential Equations", 

John Wiley and Sons, 1S64. 



[13] A. Broraan, " Introduction to Partial Differential 
Equations", Addison-Vosl sy, 1970. 

[14] F. Evans, The Physics of Fluids, Vol. 16, ho. 7, 1011, 
July 1973. 

[15] G-. Bell and S. Glasstone, "Nuclear Reactor Theory", 
Van Nostrand Reinhold Company, 1970. 

[16] R.F. Hoskins, " Generalised Functions", John Willey 
and Sons, 1979. 

[17] J.J. Duderstadt and J.L. Hamilton, "nuclear Reactor 
Physics", John Willey and Sons, 197s. 






where E can be taken as the thermal Ion energy since the 
F okker-Planck equation assumes energy independent Coulomb 
1 ogar ithm s . 



APPENDIX - II 


Properties of Dirac Delta function and Its Derivatives 

The Dirac delta function is defined as the derivative 
of the unit step function and has following properties [16,17]; 

6(x-x 0 ) = 0 , i x Q (II-1) 

= 60 , X = Xq 


6(x) = 6(-x) 


(II- 2) 


6 ( a x) = 6 (x) 

I a : 


(H-3) 


co 

/ 6 (x-x ) dx = 1 

— CO 


( II-4) 


/ f(x) 6 ( x-x Q ) dx = f(x Q ) (II-5) 

— 00 

One can differentiate the Dirac delta function as 
ui any times •' s one wishes. The u ' ' derivative has the following 
properties : 



f & W (x-a) 

CO 

f(x) dx - (-1)® 

n m i 

a f | 

, HI j 

dx 1 

x=a 

( II-6) 

6 (m) (x) = 

(-l) m 6^ m ) (-x) 


(II- 7 ) 



APPENDIX - III 


TEE UEMODBIED SECOND OIojjsH PERTURBATION EETECu 


We will start with the first order solution alven by 
Equation (3.18) in Chapter 3 . 


H) 



[6 + (l~^ 2 ) 6" I 


2p. 6 ' I] 


(III-l) 


where 


6 = 6(x q ), 


6 ' = 6(xJ , 

3 jd v O' r 


6 " £ ±— 6( Xq ) , 

3 M’ 


and other quantities are as in Chapter 3. 


The second order perturbation source will be i 

S 2= 5T i - 2,6' I] 

(III- 2) 

= ^UV)[6'+(lV)6"'l - 2H6”I 

+ 6 ,: (-2p) I - 2 6' I] 



SI 

S 2 = ~ [d-2 I) (1-j.i 2 )6 " + (4^t 3 -4n) I 6 m 

+ (1 + p 4 - 2 u 2 )I 6 1V + (1-2 I)(-2p)6 
+ (12p 2 -4)I 6 (4p 3 -4p)6 m I ] 

S I 

S 2 = TT f (1 “ 2 !) (“2J-06 T + (1 + 14 I p 2 - - P- 

+ S ( p. 3 — |_i ) I '6 + (l-!-p 4 -2p 2 ) I 6 iv ] 


So, simply Toy replacing S-^ by S 2 in Equation (3.2 
get the second order angular fZux density <f> 2 as. 

S I 

O p (x,E,p) = [ - 2p I 6 ' + 4p I 1 6 ' + (1-p 

+ (14p 2 -6) ^6*' + 8(h 3 -p)I 1 6 m 

+ (1 + p 4 -- 2p 2 )l 1 6 iv ] 


wh er e : 


I 1 = I 1 (B) 


r ' Mi ll M ill 

U(E') 


dE’ , 


- 6 1 ) 6 ' 


(I II- 2) 


ve can 


) I 6" 


(III-3) 


E 



and 1(E) is as in Eq. (3.9). 

low carrying out d|i integration, one can get the flux density 
as s 

1 

0 (x,E) = 2 tt; / (x,E,n) dp. (III-4) 

-1 

= I? [2 I - 4 I x - 2 I + 28 ^ - 43 ^ + 24 Lj 
= 0 

This agrees with the conclusion of the unmodified first order 
method, as presented in Section 3 - 1 - 
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THE LISTING OF THE COMPUTER PROGRAM 



OEU H*., *IK f MJ,T,p; f l,T,T.U 

pkm , i°)jP,TPT.r,iPt,Pi , m.Rt:, m,Ti,iPL.Tg; 

IPljTi , IPL.P1. ,1PI j T2,IPUT3,IPLR2,1PLP3 

ni"E;MSTUv ^hfioot j , a T 2 ( 1 o o 1 j 

n h-1 , 

7 = 1 . 

KK =!<■.' on 

►IS 2,4 
7,ft=J , 

MKsy.brjo 

THsrE, 

PFBiO. 

Ku=20P 
Ri P 1 s 1 

NP2sJ U 
JwM=38 . 

JJ=2 

Pi,' .,=i i ’ r ’S T K 

PJsKO+l 

YZXs^O 

Hsf PR-PMTN)/YZX 
oi=3.i*159 
F0=8 . 85'iF-l 2 
Ffc;=4.*Pi*en 
rFT=1 .860438F-27 
TYPE* , BM # Z f KS 
TYPE* ,«I,ZS,MK f TH 
TYPE* , H , FF 
ME=CFI*5.485<>7F-4. 

SE=1.6021E-19 

,0072766*CFT 
FS=ES*SpJ*l,b3 
th=th*sh;»i.f:3 
FMTN=K M JN*Sfe:*l .E3 
f*I = HT*CFT 

i,o=t./ rso«T( csf/eo) #86*47* czs+i )/th)} 

Tfc:=3.*TH*IjP*FE/(7**5K*SF) 

LlBC^^^X/fBM+^D )*ES*F£*t.O/(Z*ZS*SE*Si!:) 

Tjt*LT*TH/ES 



n = 2 . * P T * ( 8 F / F K ) * f S' F / F K 5 4 7. S * £ ? •■ Vi / -1 T ) * A T ■ 0 G C M > * S F * W F * 8 K 
r = P ],^(SF,/FF) * (Rh'/FF > *7.R*Al/lCfl,T) 
n = A*74;7, 

Rrn*7*Z 

csc*7*y. 

wsutKi^i.h:^ 

T = KO 

.A (l f K0+ !)=<).() 

A'I?lKy + nsO.O 
F = F & m t,l , 5 ♦ H 
fisP/H+JUSORTfEJ 
Fsl , /ll 
M = T + 1 

ALfDsALCMHFtH 

T=C/CE*K) 

FisT/U 

AI2(T)=AI2fM)+Fl*H 
F=F~H 
T = T-1 

TFCI.LF.O) GQTQ2G 
GOTO JO 
rONTTMOE 
XrtsAT.(t) 

TYPE*, XM 

TPIjT?sO.O 

IPM’3sO,0 

IPLR2=0 t 0 

TPLR 3*0 . 0 

TPLT 1*0,0 

XPLRlPO.O 

AIPJisQ.O 

ni'F^so.u 

AlPfaO.O 

PlPT«0.0 

Fs2 

FSTsFMT'N+H 



p IjT2-u.u 
I PT.H = U . u 
TFPRsu.O 
TPt,RT = 0,0 
TPI.HFrO.O 

TFHRTsr .0 
tfhr F rO ,0 
PHI TsO.O 
ni,j-r 1F = 0,0 
T>f,T 2 T =0 . 0 
D liT2F=f».0 
X = A L f K 1 

T •) + 0 „ 5 *H 

SO RLAsf AMT-1 ) + AM I) )*0.5 

WLA = f fr;r.-A + X)*0.5 
no ii 0 = 1 - 1 , k 

TFCWT.A.IiT, RLf J) . A*>0. WfcA „GT. AT.,( Jtl ) 1GDT025 
11 CONTINUE 

TYPE* , I , K 
7 b A0C=>,J— 1 

FlaEMlW+H*AfjC + H*f AM J)-WLA>/( AM‘ T ) -AM J + 1 ) 1 
EOlsR+M . 

KK=1~1 

FIXE?aO.O 

TFfEOl.GT.EnGOTnsi? 

F0=E+U,5*H 
«t CONTINUE, 

Ti) = C/K0**2 
UO=B/EO+A*50RT(EO) 

ELA0=AT,(KK+1) 
r>EN=2.*EI.A0-ELA 
TFa>EN t LT,0,0)TYPE*,DEA 
FIX=T0/(U0*DEN1 
pj y;;'?= ?< 'lXE2+Flx*H 

FOtsFGt+U 

KK=KKt 1 

TFfKOi .GT.EDGnTHg? 
r u =EO+H 

GOTO ft i 

v>2 CONTINUE 

Pl.f.osf f AT.fKKl + AM KK+t ) ) *0 . 5+k'li A ) *0 * 5 



? ? 0 = fa /KO+A* $ a R x CEO) 

Tfa=c/EO**2 

Pe; ,, = ?.*Er.AO-FL,.« 

T F ( l)F *i , LT . 0 . 0 ) T* Pf* f DEM 

PIXsTO/Cno*DFfO 

F] VK?=FIXE?+FIY*CE1 -E01 +H) 

T = r , /F**2 
n=n./s-+A*soRTCE) 

FT Y£=FTXF2 

PL* ISPT.11 ■+F1XE*H 

PL,T) i=p;,,tit+FXXE*{H/F)*(R/1’ ) 

dl T i. tr =PLTiF+FiyE*H* assort (F)/u 

Ai?csAT2fi-n/Aij(i-n 

TFDR=TFUR+H*A1?C 

TFORT=TFnRT + Al?C*(H/F)*C!VH) 

TFORE=TFORF+AI?C*H*A*SQRT(F)/U 

wr J A=CET,A-X)*0,S 

DO 15 JsI-l r KO 

TPf WI,A.LT.ALCJ),AHD.WLA.nT.AL(J + l ))GDT026 

1 fa CONTINUE 

’ P YPE*,T,K 
26 «KC=T-1 

Fl=E’UN + H*ABCtH*C A!.( J)-NT.A)/f AT.,( J) -AFj ( J+t ) ) 

FlXE?eO.O 

FOlsFS-H 

KKsKO+1 

TF(E01,faT,Fl)GOTD84 
F0=ES-O,5*H 
«3 T0=C/E0**2 

UO=B/E0+A*SQRT(E0) 

E(iA0=0.5*C»L(KK)+Ar.(KK*n) 

DFW=t?.*E?,AO-FL» 

DEM = -i)FiX 

TFtDFH.LT.O.O)TYPE* f OEN 

FIXxTO/CUOIDFN) 

FIXE?=FIXF2+FIX*H 

F01=F01 -U 
KK=KK-1 

TFfEOJ ,LT.Fl)GOT r >84 



T**! i) ss ^, r f! m 
fTOTuQ 3 

*fi HjVTTftf T f. 

Ffj a v=r C a t. C K Jf } + *r,A)*o, *? 

FO=(FOl + H + F1 1 *0,5 
no = A/!?:O+B.*R 0 iPXCEO} 

tu=c/ep **2 

1 nt;»is?, *j£TiAO-KL,A 
r >F M :=’»UF^ 

T F ( UFA, , |jT . 0 » 0 ) ^ V’Pf;* , 

TF f PFU , l.F ,0.0) CJQTf.jR 5 
F I Vs: TO / ( M 0 * r> F W ) 

9 5 Frx= 0„0 

FIXEt?ssFIVK? + Fiy*fE')i+H-E1 ) ■ 

fixe;=ftxf2 

PLT2 = PT,I2+FIXe:*H 

PLT2I = PLT2I + FIXe*{H/F)*(P/T!) 

PLT2F=PL r2F+FIXE»H* ASSORT (P)/Ii 
»ljTAst./AL<I-l) 

TPt.R = IPkR + 0.5*H*XH*AT,IA/FS 
TPT,RF3TPT,RE+0.5*H*XM*C AI j TA/ESJ*A*SOHT(F)/U 
TPT.«TsTPTiRT + 0.5MH/£)*X<l*C AI.,TA/ES)*CB/»>) 

FsF + H 

TFfI.EO,K)CnT060 

TsT + 1 

OQT050 

fiU CONTIMUF 

X l! = X / X M 
X=ftL(K.-X) 

mrrvx ( »^c*np + pr \ 1+PLT2 ) *Xrt/RS 
CHFCTs* C«TFORT+Pl<T 1 T+PLI2T)*X M /F5 
ChFCF" ( “' T ’F n RF+PLI 1 F+PLI2F) *XM,/FS 
Tpr.TsTPi.n+rH^CK 
TxTPTjRT +CHFCI 
TpT.X”- T P r -RF4 ri ’FCF 
ALT-DsAt-CK-t )-Al.(F) 
PlM=2.*X*tTpr,T1-TPT J T)/A|jT J D 
»lMl = 2.*V*f 1PLT1-IPI.T2VAH.D 
Fi"i:=2.*X*(lPl,T’f>XP3UT3)/Ar J T,D 
PX"J T=2.*X*CTPLRr-lPr J R2)/Ar J l,n 
t»T W t 5-s? - *X* f TPLRF-Tpr.H 3 ) / hhW 



p t^U“2.*X*aPLRl«lPI ) R)/AU g n 

PX M ls5-PI«t 
TPT/J’I sTPTjT 
Tpr,T?=Tpr,:cT 
TPVnsTPM’P- 
TPLR2=TPT.«T 
Tpr.R3sTPl.RE 
TPLRIsTPTjR 
ESK=:EST/(SK*1 .Eil 

T F f K , Pin , J J ) TYPE* , K , X” , I PT.R , 1 PijP l , 1 Pt.PfS 

TFfK.fcO.JJ VrYPJT*, K,fcSK, IP1,T, Tpr.'l'T , TPT.TE 

T r ( K . R 0 „ ,T j ■) T Y p K * , P T f n , P J m 1 T , P I M 1 E 

IF ( K . EO. JJ ? TYPE* , PTW , PT .'4 T , Pi NK 

TFfK,EO,JJM,.i=.lj + PPl 

TFCK.EO, s»P2 

A lPLsAlPL+TPT.«*A!jLn/XM 

PrPL = BTpI, + iPTiT*Ar>Ln/XM 

A % P t s A T p ■ T + p I N J * A ? , I , ! V X H 

^ X p X = B T p T + P x w * A 1 , r , D / X M 

KssK + 1 

FSTsEST+H 

TFCK.LE.KJ) GDT07 0 

TYPE*, AIPL, BIPli, AXPT,BJPT 

STOP 

END 



A 70565 




